GÉO-VISUALISATION AVEC R

Action nationale de formation
MATE SHS

Cette formation de 3 heures 30 porte sur la visualisation de données géographiques sous R. Y sont abordées les traitements SIG de base, la cartographie thématique (cartes en figurés proportionnels, cartes choroplèthes, cartes de typologie, etc) et des cartes reposant sur des techniques plus avancées comme les cartes sur grille ou les cartes de discontinuité.

Sont ici principalement utilisés les packages sf (manipulation de données spatiales) et cartography (cartographie thématique). Chaque exemple propose une chaine de traitement depuis le chargement des données, leur mise en forme, jusqu’à la mobilisation de méthodes adaptées permettant de répondre à des questions spatiales.

Les exemples proposés sont basés sur différentes échelles, du local (région Occitanie) en passant par l’Europe et le Monde.

Pré-réquis

Les packages

Sont décrits ici la version de R avec laquelle a été produit ce support de formation ainsi que les packages qui seront utilisés dans la formation et leur version respective.

library(devtools)
session_info(pkgs = c("cartography","sf"), include_base = FALSE)
## Session info -------------------------------------------------------------
##  setting  value                       
##  version  R version 3.4.4 (2018-03-15)
##  system   x86_64, linux-gnu           
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  tz       Europe/Paris                
##  date     2018-09-27
## Packages -----------------------------------------------------------------
##  package     * version date       source        
##  abind         1.4-5   2016-07-21 CRAN (R 3.4.4)
##  cartography   2.1.2   2018-09-18 CRAN (R 3.4.4)
##  class         7.3-14  2015-08-30 CRAN (R 3.4.2)
##  classInt      0.2-3   2018-04-16 CRAN (R 3.4.4)
##  curl          3.2     2018-03-28 CRAN (R 3.4.4)
##  DBI           1.0.0   2018-05-02 CRAN (R 3.4.4)
##  digest        0.6.17  2018-09-12 CRAN (R 3.4.4)
##  e1071         1.7-0   2018-07-28 CRAN (R 3.4.4)
##  graphics    * 3.4.4   2018-04-21 local         
##  grDevices   * 3.4.4   2018-04-21 local         
##  grid          3.4.4   2018-04-21 local         
##  httr          1.3.1   2017-08-20 CRAN (R 3.4.4)
##  jpeg          0.1-8   2014-01-23 CRAN (R 3.4.4)
##  jsonlite      1.5     2017-06-01 CRAN (R 3.4.4)
##  lattice       0.20-35 2017-03-25 CRAN (R 3.4.2)
##  magrittr      1.5     2014-11-22 CRAN (R 3.4.4)
##  MASS          7.3-49  2018-02-23 CRAN (R 3.4.3)
##  methods     * 3.4.4   2018-04-21 local         
##  mime          0.5     2016-07-07 CRAN (R 3.4.4)
##  openssl       1.0.2   2018-07-30 CRAN (R 3.4.4)
##  plyr          1.8.4   2016-06-08 CRAN (R 3.4.4)
##  png           0.1-7   2013-12-03 CRAN (R 3.4.4)
##  prettymapr    0.2.2   2017-09-20 CRAN (R 3.4.4)
##  R6            2.2.2   2017-06-17 CRAN (R 3.4.4)
##  raster        2.6-7   2017-11-13 CRAN (R 3.4.4)
##  Rcpp          0.12.18 2018-07-23 CRAN (R 3.4.4)
##  rgdal         1.3-4   2018-08-03 CRAN (R 3.4.4)
##  rgeos         0.3-28  2018-06-08 CRAN (R 3.4.4)
##  rjson         0.2.20  2018-06-08 CRAN (R 3.4.4)
##  rosm          0.2.2   2017-04-07 CRAN (R 3.4.4)
##  sf            0.6-3   2018-05-17 CRAN (R 3.4.4)
##  sp            1.3-1   2018-06-05 CRAN (R 3.4.4)
##  spData        0.2.9.4 2018-09-15 CRAN (R 3.4.4)
##  stats       * 3.4.4   2018-04-21 local         
##  tools         3.4.4   2018-04-21 local         
##  units         0.6-1   2018-09-21 CRAN (R 3.4.4)
##  utils       * 3.4.4   2018-04-21 local

Pour exécuter les programmes proposés par la formation, installez / mettez à jour R (version 3.5.1 minimum), R Studio. Lancez ensuite R studio et installez les packages au moyen des commandes suivantes (prévoir bien 10 minutes)

# install.packages("sf")
# install.packages("devtools")
# install.packages("cartography")
# install.packages("countrycode")
# install.packages("readxl")
# install.packages("animation")
# install.packages("rmapshaper")
# install.packages("rnaturalearthdata")
# install.packages("rnaturalearth")
# install.packages("cartogram")
# install.packages("reshape2")
# install.packages("eurostat")
  • Le package sf est un package qui permet d’importer, gérer et transformer des données géographiques (gestion des projections, opérations SIG).
  • Le package cartography permet de créer et intégrer des cartes thématiques dans sa chaine de traitement en R. Il permet des représentations cartographiques tels que les cartes en symboles proportionnels, des cartes choroplèthes, des typologies, des cartes de flux ou des cartes de discontinuité. Il offre également des fonctions qui permettent d’améliorer la réalisation de la carte, comme des palettes de couleur, des éléments d’habillage (échelle, flèche du nord, titre, légende…), d’y rattacher des labels ou d’accéder à des APIs cartographiques.
  • Le package…

Téléchargez les données

Téléchargez les données sur la page github. Décompressez les données sur votre espace de travail.

Lien github à ajouter

Exo 1 - Prise en main

Commandes de base

Définissez votre répertoire de travail.

setwd("Votre_repertoire_de_travail")

Consulter le contenu. Et regarder ce qu’il y a dans le répertoire “data”

list.files()
## [1] "data"        "data.zip"    "index_files" "index.html"  "index.Rmd"  
## [6] "outputs"     "pics"        "prg"
list.files("data")
## [1] "download"  "Europe"    "France"    "Hérault"   "Occitanie" "world"

Import d’une couche SIG dans R avec le package sf. Pacakge basé sur GEOS (Geometry Engine Open Source) et GDAL (Geospatial Data Abstraction Library).

library("sf")
## Linking to GEOS 3.6.2, GDAL 2.2.3, proj.4 4.9.3
communes <- st_read(dsn = "data/Hérault/communes.shp", stringsAsFactors = F)
## Reading layer `communes' from data source `/home/nlambert/ownCloud/ANF2018 - R geoviz/data/Hérault/communes.shp' using driver `ESRI Shapefile'
## Simple feature collection with 343 features and 5 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 662675 ymin: 6234874 xmax: 796333 ymax: 6319348
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs

Voir la table atributaire

communes
## Simple feature collection with 343 features and 5 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 662675 ymin: 6234874 xmax: 796333 ymax: 6319348
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
## First 10 features:
##    INSEE_COM                       NOM_COMM           STATUT SUPERFICIE
## 1      34029                        BELARGA   Commune simple        415
## 2      34290 SAINT-VINCENT-DE-BARBEYRARGUES   Commune simple        222
## 3      34032                        BEZIERS   Sous-prfecture       9567
## 4      34082                    COMBAILLAUX   Commune simple        902
## 5      34108                     FRONTIGNAN Chef-lieu canton       3984
## 6      34166                      MONTBLANC   Commune simple       2716
## 7      34066                    CAZEVIEILLE   Commune simple       1619
## 8      34296                      SAUSSINES   Commune simple        630
## 9      34269        SAINT-JEAN-DE-MINERVOIS   Commune simple       3277
## 10     34012                     ARGELLIERS   Commune simple       5065
##    NOM_DEPT                       geometry
## 1   HERAULT POLYGON ((742075 6272005, 7...
## 2   HERAULT POLYGON ((770896 6289405, 7...
## 3   HERAULT POLYGON ((718759 6244261, 7...
## 4   HERAULT POLYGON ((761606 6284489, 7...
## 5   HERAULT POLYGON ((758738 6257679, 7...
## 6   HERAULT POLYGON ((728067 6248191, 7...
## 7   HERAULT POLYGON ((762395 6294296, 7...
## 8   HERAULT POLYGON ((784582 6294234, 7...
## 9   HERAULT POLYGON ((686897 6252592, 6...
## 10  HERAULT POLYGON ((756922 6287661, 7...
head(communes)
## Simple feature collection with 6 features and 5 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 710431 ymin: 6244261 xmax: 771300 ymax: 6292043
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
##   INSEE_COM                       NOM_COMM           STATUT SUPERFICIE
## 1     34029                        BELARGA   Commune simple        415
## 2     34290 SAINT-VINCENT-DE-BARBEYRARGUES   Commune simple        222
## 3     34032                        BEZIERS   Sous-prfecture       9567
## 4     34082                    COMBAILLAUX   Commune simple        902
## 5     34108                     FRONTIGNAN Chef-lieu canton       3984
## 6     34166                      MONTBLANC   Commune simple       2716
##   NOM_DEPT                       geometry
## 1  HERAULT POLYGON ((742075 6272005, 7...
## 2  HERAULT POLYGON ((770896 6289405, 7...
## 3  HERAULT POLYGON ((718759 6244261, 7...
## 4  HERAULT POLYGON ((761606 6284489, 7...
## 5  HERAULT POLYGON ((758738 6257679, 7...
## 6  HERAULT POLYGON ((728067 6248191, 7...
View(communes)

L’instruction plot permet d’afficher la couche. L’instruction st_geometry permet d’accéder à la variable définissant les géométries.

plot(st_geometry(communes))

Afficher la couche avec des parametres graphiques

plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)

Introduction au SIG avec R / Manipuler les informations spatiales

Nous cherchons à identifier les communes qui ont une partie de leur territoire situé à moins de 20 km du centre de Sete.

Nous commencçons par extraire la commune de Sète.

macommune <- "SETE"
monpoly <- communes[communes$NOM_COMM == macommune,]
plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)
plot(st_geometry(monpoly), col="red", border="purple", lwd=1, add=T)

Extraitre le centroide de la commune de Sete.

moncentre <- st_centroid(x = monpoly)
## Warning in st_centroid.sf(x = monpoly): st_centroid assumes attributes are
## constant over geometries of x
plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)
plot(st_geometry(monpoly), col="red", border="purple", lwd=1, add=T)
plot(st_geometry(moncentre), pch=20, col="black", cex=2, add=T)

Calculer une zone tampon

mydist <- 20000

buff <- st_buffer(x = st_geometry(moncentre), dist=mydist)
plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)
plot(st_geometry(monpoly), col="red", border="purple", lwd=1, add=T)
plot(st_geometry(moncentre), pch=20, col="black", cex=2, add=T)
plot(buff, col=NA, border="black", lty=2,add=T)

Récupérer la liste des communes

communes$buff <- st_intersects(st_geometry(communes), st_geometry(buff), sparse = FALSE)
head(communes)
## Simple feature collection with 6 features and 6 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 710431 ymin: 6244261 xmax: 771300 ymax: 6292043
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
##   INSEE_COM                       NOM_COMM           STATUT SUPERFICIE
## 1     34029                        BELARGA   Commune simple        415
## 2     34290 SAINT-VINCENT-DE-BARBEYRARGUES   Commune simple        222
## 3     34032                        BEZIERS   Sous-prfecture       9567
## 4     34082                    COMBAILLAUX   Commune simple        902
## 5     34108                     FRONTIGNAN Chef-lieu canton       3984
## 6     34166                      MONTBLANC   Commune simple       2716
##   NOM_DEPT                       geometry  buff
## 1  HERAULT POLYGON ((742075 6272005, 7...  TRUE
## 2  HERAULT POLYGON ((770896 6289405, 7... FALSE
## 3  HERAULT POLYGON ((718759 6244261, 7... FALSE
## 4  HERAULT POLYGON ((761606 6284489, 7... FALSE
## 5  HERAULT POLYGON ((758738 6257679, 7...  TRUE
## 6  HERAULT POLYGON ((728067 6248191, 7... FALSE

Récupérer le nombre de communes à moins de 20 km

nb <- dim(communes[communes$buff == T,])
nb[1]
## [1] 40

Extraction et affichage des communes

communes20km <- communes[communes$buff == TRUE,]
plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)
plot(st_geometry(monpoly), col="red", border="purple", lwd=1, add=T)
plot(st_geometry(moncentre), pch=20, col="black", cex=2, add=T)
plot(buff, col=NA, border="black", lty=2,add=T)
plot(st_geometry(communes20km), col=NA, lwd=2, border="red", add=T)

Ajout d’un titre

plot(st_geometry(communes), col="#aec8f2", border="darkblue", lwd=1)
plot(st_geometry(monpoly), col="red", border="purple", lwd=1, add=T)
plot(st_geometry(moncentre), pch=20, col="black", cex=2, add=T)
plot(buff, col=NA, border="black", lty=2,add=T)
plot(st_geometry(communes20km), col=NA, lwd=2, border="red", add=T)
montitre <- paste0("Il y a ", nb[1], " communes situées à moins de ", mydist/1000, "km de ",tolower(macommune))
title(montitre)

Export (si besoin)

st_write(obj = communes20km, dsn = "outputs/resultat.shp")

A vous de jouer !

Reproduisez cette procédure avec la commune et la distance de votre choix.

Pour accéder au listing des communes d’Hérault

communes$NOM_COMM

Exo 2 - Cartes en figurés proportionnels

Objectifs :

  • Thématique : Réaliser une carte des catégories Socio professionnelles.
  • Espace d’étude : Occitanie, communes
  • Visualisation associée : carte en figurés proportionnels (+ cartogrammes de Dorling)
  • packages utilisés : sf, cartography, cartogram, readxl

Créer un nouveau fichier “Exo2.R”.

Chargement et préparation des données

Pour bien commencer, il est fondamental de jeter un oeil au fichier d’entrée pour comprendre son organisation et savoir comment l’importer sous R.

Ouvrez le fichier data/France/INSEE/pop-act2554-csp-cd-6814.xls Nous nous intéressons ici à l’onglet “COM_2014”.

library("sf")
library("readxl")

Visualisation : carte en figuré proportionnel

library("cartography")

variable visuelle taille adaptée aux caractères quantitatifs absolus.

Visualisation : carte combinant figuré proportionnel et typologie

variable visuelle taille + variable visuelle couleur

Visualisation : cartogramme de Dorling

library("cartogram")

Exo 3 -

Chargement des packages nécessaires

xxx

xxx

Exo 4 - Cartes de discontinuité

Objectifs :

  • Thématique : Les inégalités d’IDH et/ou d’espérance de vie
  • Espace d’étude : Le Monde, pays.
  • Visualisation associée : carte chorplèthe, discontinuités
  • packages utilisés : sf, cartography, rnaturalearth, readxl, countrycode

Import du fond de carte

méthode 1 : à la main (après téléchargement sur le site Natural Earth)

library("sf")
countries <- st_read(dsn = "data/world/naturalearth/ne_110m_admin_0_countries.shp", stringsAsFactors = F)  

Méthode 2 : directement via R

library("sf")
url <- "https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/cultural/ne_110m_admin_0_countries.zip"
download.file(url = url, destfile = "data/download/countries.zip")
unzip("data/download/countries.zip", exdir = "data/download", unzip = "internal")
file.remove("data/download/countries.zip")
countries <- st_read(dsn = "data/download/ne_110m_admin_0_countries.shp", stringsAsFactors = F)

Méthode 3 : directement via R

library("rnaturalearth")
countries <- ne_countries(scale = 50, type = "countries", continent = NULL,
                          country = NULL, geounit = NULL, sovereignty = NULL,
                          returnclass = "sf")

Le fond de carte peut être visualisé avec l’instruction plot

plot(st_geometry(countries))

On ne conserve que les champs utiles. On les renomme.

countries <- countries[,c("adm0_a3", "admin", "geometry")]
colnames(countries) <- c("id","name","geometry")

On Supprime l’Antarctique

countries <- countries[countries$id != "ATA",]

On supprime les lignes vides

countries <- countries[!is.na(countries$id),]
plot(st_geometry(countries))

Import des données atributaires

Import des données sur l’éspérance de vie (source : Banque mondiale)

library("readxl")
file <-"data/world/worldbank/API_SP.DYN.LE00.IN_DS2_en_excel_v2_10081535.xls"
sheet <- "Data"
lifexp <- data.frame(read_excel(file, skip = 3, sheet = sheet))
lifexp <- lifexp[,c("Country.Code","Country.Name","X2016")]
colnames(lifexp) <- c("id","name","lifexp2016")
lifexp$lifexp2016 <- as.numeric(as.character(lifexp$lifexp2016))

Import des données sur l’Indice de deveoppement humain (source : Human Development Report)

file <-"data/world/hdr/Human Development Index (HDI).csv"
hdi <- read.csv2(file = file, sep = ",",skip = 1)
hdi <- hdi[,c("Country","X2016")]
head(hdi,4)
##        Country X2016
## 1  Afghanistan 0.494
## 2      Albania 0.782
## 3      Algeria 0.753
## 4      Andorra 0.856

Le tableau de données sur l’IDH ne possède pas de codes géographiques, mais uniquement le nom des pays. Cela rend la jointure avec le fond de carte compliquée. Une solution : le package countrycode

library("countrycode")
hdi$id <- countrycode(hdi$Country, 'country.name', 'iso3c')
## Warning in countrycode(hdi$Country, "country.name", "iso3c"): Some values were not matched unambiguously:  Eswatini (Kingdom of)

Un peu de mise en forme. Les codes sont bien là. C’est presque magique… Mais attention, dans la pratique, une vérification manuelle s’impose.

hdi <- hdi[,c("id","Country","X2016")]
colnames(hdi) <- c("id","name","hdi2016")
hdi$hdi2016 <- as.numeric(as.character(hdi$hdi2016))  
head(hdi,4)
##    id         name hdi2016
## 1 AFG  Afghanistan   0.494
## 2 ALB      Albania   0.782
## 3 DZA      Algeria   0.753
## 4 AND      Andorra   0.856

Cartographie

Chargement du package

library("cartography")

Changer la projection du fond de carte (voir le site http://spatialreference.org)

countries <- st_transform(countries, "+proj=cea +lon_0=0 +lat_ts=45 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs ")
plot(st_geometry(countries))

Jointures des deux fichiers de données

countries <- merge(countries, lifexp, by="id", all.x=TRUE)
countries <- merge(countries, hdi, by="id", all.x=TRUE)
countries <- countries[,c("id","name.x","lifexp2016","hdi2016","geometry")]
colnames(countries)[2] <- "name"

On travaille sur l’indice de développement humain

analyse de la distribution

var <- countries$hdi2016
hist(countries$hdi2016, probability = TRUE, nclass = 30)

On choisit une discretisation par la méthode des quantiles sans classe centrale

breaks <- getBreaks(v = var, nclass = 6, method = "quantile")
cols <- carto.pal(pal1="blue.pal", n1 = 3, pal2 = "orange.pal", n2 = 3, middle = FALSE, transparency = FALSE)

Réalisation de la carte

layoutLayer(extent = countries,
            bg="#EEEEEE",
            title = "Indice de développement humain en 2016",
            scale = NULL,
            sources = "Human development Report, 2018",
            author = "les serial mappers"
            )
choroLayer(x = countries, 
           var = "hdi2016",
           breaks = breaks,
           col = cols,
           border = "white",
           lwd=0.2,
           add = TRUE,
           legend.pos = "topleft",
           legend.title.txt = "IDH, 2016",
           legend.values.rnd = 2)

Et si on ajoutait des lignes de discontinuités à cette carte?

L’instruction getBorders (du package cartography) permet d’extraire les frontières entre les unités spatiales.

borders <- getBorders(countries)
plot(st_geometry(borders), col="black", lwd=1)

Quid des discontinuités paritimes ? Par exemple entre l’Italie et la Tunisie ? Ou entre le Maroc et l’Espagne ? L’instruction getOuterBorders (package cartography) permet de les générer. Cette instruction peut prendre un peu de temps.

outer <- getOuterBorders(x = countries, res = 100000, width=500000)
plot(st_geometry(borders), col="#CCCCCC", lwd=1)
plot(st_geometry(outer), col="red", lwd=1, add=T)

On assemble les deux couches avec rbind

b <- rbind(borders,outer)

Tout est prêt pour réaliser une carte de discontinuités

# 1 - la carte choroplèthe
layoutLayer(extent = countries,
            bg="#EEEEEE",
            title = "Indice de développement humain en 2016",
            scale = NULL,
            sources = "Human development Report, 2018",
            author = "les serial mappers"
            )
choroLayer(x = countries, 
           var = "hdi2016",
           breaks = breaks,
           col = cols,
           border = "white",
           lwd=0.2,
           add = TRUE,
           legend.pos = "topleft",
           legend.title.txt = "IDH, 2016",
           legend.values.rnd = 2)

# 2 - les discontinuités

discLayer(x = b, df = countries,
          var = "hdi2016", col="black", nclass=4,
          method="quantile", threshold = 0.25, sizemin = 0.2,
          sizemax = 8, type = "abs", 
          legend.title.txt = "Discontinuities\nabsolues",
          legend.pos = "bottomleft", add=TRUE)

BONUS

Utiliser des données OSM

Ressources

  • BEGUIN Michèle, PUMAIN Denise, La représentation des données géographiques, Statistique et cartographie, coll. Cursus, Armand Colin, nouvelle édition 2000, 192p.
  • BERTIN Jacques, Sémiologie graphique, Monton-Gauthier-Villars, 1967, 1973, 432 p. Disponible à ce jour dans la collection Réimpression de l’EHESS, 1998.
  • BERTIN Jacques, La graphique et le traitement graphique de l’information, Flammarion, 1977, 280 p.
  • BRUNET Roger, La carte mode d’emploi, Fayard-Reclus, 1987, 270p.
  • ElementR (Groupe), R et espace, traitement de l’information géographique, 2014, Framabook.
  • LAMBERT Nicolas : Carnet de recherche neocarto : https://neocarto.hypotheses.org
  • LAMBERT Nicolas, ZANIN Christine : « Manuel de cartographie. Principes, méthodes, applications », cursus, Armand Colin, 18 mai 2016, 224p.
  • Giraud Timothée : Carnet de recherche rgeomatic : https://rgeomatic.hypotheses.org/
  • ZANIN Christine, TREMELO Marie-Laure, Savoir faire une carte, Aide à la conception et à la réalisation d’une carte thématique univariée, coll. Sup géographie, Belin, 2003, 200p.

Nicolas Lambert, Ronan Ysebaert

Sète, 16 nov. 2018